home <- here::here()
library(tidyverse)
library(tidylog)
library(RCurl)
library(sdmTMB)
library(RColorBrewer)
library(devtools)
library(patchwork)
library(ggstats)
library(ggh4x)
library(sdmTMBextra)
# Source map-plot
#source_url("https://raw.githubusercontent.com/maxlindmark/cod-interactions/main/R/functions/map-plot.R")
source(paste0(home, "/R/functions/map-plot.R"))Food competition analysis
Load packages
Read data & prepare data
d <- read_csv(paste0(home, "/data/clean/aggregated_stomach_data.csv")) %>%
drop_na(group) %>%
drop_na(oxy)
t <- d %>% filter(fle_kg_km2 > 0)
t <- d %>% filter(mcod_kg_km2 > 0)
t <- d %>% filter(scod_kg_km2 > 0)
# Calculate relative prey weights (saduria and benthos)
d <- d %>%
rename(oxygen = oxy) %>%
mutate(tot_weight = rowSums(select(., ends_with('_tot'))),
benthic_weight = amphipoda_tot + bivalvia_tot + gadus_morhua_tot +
gobiidae_tot + mysidae_tot + non_bio_tot +
other_crustacea_tot + other_tot + other_pisces_tot + platichthys_flesus_tot +
polychaeta_tot + saduria_entomon_tot) %>%
rename(saduria_weight = saduria_entomon_tot,
flounder_density = fle_kg_km2,
large_cod_density = mcod_kg_km2,
small_cod_density = scod_kg_km2) %>%
mutate(tot_rel_weight = tot_weight / (pred_weight_g - tot_weight),
benthic_rel_weight = benthic_weight / (pred_weight_g - tot_weight),
saduria_rel_weight = saduria_weight / (pred_weight_g - tot_weight)) %>%
dplyr::select(-ends_with("_tot")) %>%
dplyr::select(-predator_latin_name, date) %>%
# Add small constant to large cod density because we want to take the log of it
mutate(large_cod_density = ifelse(large_cod_density == 0,
min(filter(d, mcod_kg_km2 > 0)$mcod_kg_km2)*0.5,
large_cod_density),
flounder_density = ifelse(flounder_density == 0,
min(filter(d, fle_kg_km2 > 0)$fle_kg_km2)*0.5,
flounder_density),
small_cod_density = ifelse(small_cod_density == 0,
min(filter(d, scod_kg_km2 > 0)$scod_kg_km2)*0.5,
small_cod_density)) %>%
# scale variables
mutate(fyear = as.factor(year),
fquarter = as.factor(quarter),
fhaul_id = as.factor(haul_id),
depth_sc = as.numeric(scale(depth)),
oxygen_sc = as.numeric(scale(oxygen)),
density_saduria_sc = as.numeric(scale(density_saduria)),
flounder_density_sc = as.numeric(scale(log(flounder_density))),
large_cod_density_sc = as.numeric(scale(log(large_cod_density))),
small_cod_density_sc = as.numeric(scale(log(small_cod_density))))Quick explore
Correlation between variables
# Plot correlation between variables
d_cor <- d %>%
dplyr::select("oxygen_sc", "density_saduria_sc", "flounder_density_sc",
"large_cod_density_sc", "small_cod_density_sc", "depth_sc")
corr <- round(cor(d_cor), 1)
library(ggcorrplot)
ggcorrplot(corr, hc.order = TRUE, type = "lower", lab = TRUE, lab_size = 2.5) +
theme_classic() +
theme(axis.text.x = element_text(angle = 90, vjust = 0.3))
Sample size
d %>%
group_by(species) %>%
summarise(n = n())group_by: one grouping variable (species)summarise: now 2 rows and 2 columns, ungrouped# A tibble: 2 × 2
species n
<chr> <int>
1 Cod 5103
2 Flounder 3851d %>%
group_by(species, quarter) %>%
summarise(n = n())group_by: 2 grouping variables (species, quarter)summarise: now 4 rows and 3 columns, one group variable remaining (species)# A tibble: 4 × 3
# Groups: species [2]
species quarter n
<chr> <dbl> <int>
1 Cod 1 2906
2 Cod 4 2197
3 Flounder 1 2081
4 Flounder 4 1770Fit models
Groups are: small cod, large cod and flounder. Response variables are: saduria_rel_weight, benthic_rel_weight or total weight. The latter is only for adult cod, because essentially all prey are benthic for small cod and flounder.
Model random effect structure is selected with AIC (see script 02-)
# This is the reason we don't do total weight for flounder and small cod
d %>%
filter(tot_rel_weight > 0) %>%
group_by(group) %>%
mutate(ben_prop = benthic_rel_weight / tot_rel_weight) %>%
summarise(mean_ben_prop = mean(ben_prop))# A tibble: 3 × 2
group mean_ben_prop
<chr> <dbl>
1 flounder 0.978
2 large cod 0.592
3 small cod 0.956Covariates are: ~ 0 + fyear + fquarter + depth_sc + spatial + spatiotemporal random fields + density covariates. For saduria, we use saduria also in interaction with cod and flounder. For cod we use small cod because large and small cod are very correlated. For benthic and total prey, we instead use oxygen, more as a proxy, as the interaction variable
Main models
pred_flounder_sad <- list()
pred_flounder_ben <- list()
pred_cod_sad <- list()
pred_cod_ben <- list()
coef_sad <- list()
coef_ben <- list()
res_sad <- list()
res_ben <- list()
random_sad <- list()
random_ben <- list()
range_sad <- list()
range_ben <- list()
for(i in unique(d$group)) {
dd <- filter(d, group == i)
mesh <- make_mesh(dd,
xy_cols = c("X", "Y"),
cutoff = 5)
ggplot() +
inlabru::gg(mesh$mesh) +
coord_fixed() +
geom_point(aes(X, Y), data = dd, alpha = 0.2, size = 0.5) +
annotate("text", -Inf, Inf, label = paste("n knots =", mesh$mesh$n), hjust = -0.3, vjust = 3) +
labs(x = "Easting (km)", y = "Northing (km)")
ggsave(paste0(home, "/figures/supp/mesh_", i, ".pdf"), width = 17, height = 17, units = "cm")
# Saduria model
m_sad <- sdmTMB(saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc +
small_cod_density_sc*density_saduria_sc +
flounder_density_sc*density_saduria_sc,
data = dd,
mesh = mesh,
family = tweedie(),
spatiotemporal = "IID",
spatial = "on",
time = "year")
print(i)
sanity(m_sad)
print(m_sad)
# Benthic model
if(unique(dd$group) %in% c("large cod", "small cod")) {
m_ben <- sdmTMB(benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc +
small_cod_density_sc*oxygen_sc + flounder_density_sc*oxygen_sc,
data = dd,
mesh = mesh,
family = tweedie(),
spatiotemporal = "IID",
spatial = "off",
time = "year")
print(i)
sanity(m_ben)
print(m_ben)
} else {
m_ben <- sdmTMB(benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc +
small_cod_density_sc*oxygen_sc + flounder_density_sc*oxygen_sc,
data = dd,
mesh = mesh,
family = tweedie(),
spatiotemporal = "IID",
spatial = "on",
time = "year")
print(i)
sanity(m_ben)
print(m_ben)
}
# Spatial and spatiotemporal random effects
d_haul <- dd %>%
distinct(haul_id, .keep_all = TRUE)
preds_sad <- predict(m_sad, newdata = d_haul)
preds_ben <- predict(m_ben, newdata = d_haul)
random_sad[[i]] <- preds_sad
random_ben[[i]] <- preds_ben
# Residuals
samps <- sdmTMBextra::predict_mle_mcmc(m_sad, mcmc_iter = 401, mcmc_warmup = 400)
mcmc_res <- residuals(m_sad, type = "mle-mcmc", mcmc_samples = samps)
dd$res <- as.vector(mcmc_res)
res_sad[[i]] <- dd
samps <- sdmTMBextra::predict_mle_mcmc(m_ben, mcmc_iter = 401, mcmc_warmup = 400)
mcmc_res <- residuals(m_ben, type = "mle-mcmc", mcmc_samples = samps)
dd$res <- as.vector(mcmc_res)
res_ben[[i]] <- dd
# Ranges
range_sad[[i]] <- tidy(m_sad, effects = "ran_pars") %>% filter(term == "range") %>% mutate(group = i, model = "saduria")
range_ben[[i]] <- tidy(m_ben, effects = "ran_pars") %>% filter(term == "range") %>% mutate(group = i, model = "benthos")
# Conditional effects: flounder
nd_flounder <- data.frame(expand_grid(
density_saduria_sc = c(quantile(d$density_saduria_sc, probs = 0.05),
mean(d$density_saduria_sc),
quantile(d$density_saduria_sc, probs = 0.95)),
flounder_density_sc = seq(quantile(dd$flounder_density_sc, probs = 0.05),
quantile(dd$flounder_density_sc, probs = 0.95),
length.out = 50))) %>%
mutate(year = 2020,
fyear = as.factor(2020),
fquarter = as.factor(1),
oxygen_sc = 0,
depth_sc = 0,
small_cod_density_sc = 0)
preds_flounder_sad <- predict(m_sad, newdata = nd_flounder, re_form = NA, re_form_iid = NA, se_fit = TRUE)
preds_flounder_ben <- predict(m_ben, newdata = nd_flounder, re_form = NA, re_form_iid = NA, se_fit = TRUE)
pred_flounder_sad[[i]] <- preds_flounder_sad %>% mutate(group = i, xvar = "flounder")
pred_flounder_ben[[i]] <- preds_flounder_ben %>% mutate(group = i, xvar = "flounder")
# Conditional effects: cod
nd_cod <- data.frame(expand_grid(
density_saduria_sc = c(quantile(d$density_saduria_sc, probs = 0.05),
quantile(d$density_saduria_sc, probs = 0.95)),
small_cod_density_sc = seq(quantile(dd$small_cod_density_sc, probs = 0.05),
quantile(dd$small_cod_density_sc, probs = 0.95),
length.out = 50))) %>%
mutate(year = 2020,
fyear = as.factor(2020),
fquarter = as.factor(1),
oxygen_sc = 0,
depth_sc = 0,
flounder_density_sc = 0) #
preds_cod_sad <- predict(m_sad, newdata = nd_cod, re_form = NA, re_form_iid = NA, se_fit = TRUE)
preds_cod_ben <- predict(m_ben, newdata = nd_cod, re_form = NA, re_form_iid = NA, se_fit = TRUE)
pred_cod_sad[[i]] <- preds_cod_sad %>% mutate(group = i, xvar = "cod")
pred_cod_ben[[i]] <- preds_cod_ben %>% mutate(group = i, xvar = "cod")
# Coefficients
coefs_sad <- bind_rows(tidy(m_sad, effects = "fixed", conf.int = TRUE)) %>%
mutate(species = "Cod (m)",
response = "Saduria",
sig = ifelse(estimate > 0 & conf.low > 0, "Y", "N"),
sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig))
coefs_ben <- bind_rows(tidy(m_ben, effects = "fixed", conf.int = TRUE)) %>%
mutate(species = "Cod (m)",
response = "Saduria",
sig = ifelse(estimate > 0 & conf.low > 0, "Y", "N"),
sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig))
coef_sad[[i]] <- coefs_sad %>% mutate(group = i)
coef_ben[[i]] <- coefs_ben %>% mutate(group = i)
}filter: removed 5,698 rows (64%), 3,256 rows remainingLoading required namespace: INLA[1] "large cod"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably largeSpatiotemporal model fit by ML ['sdmTMB']
Formula: saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc +
Formula: small_cod_density_sc * density_saduria_sc + flounder_density_sc *
Formula: density_saduria_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
coef.est coef.se
fyear2015 -8.63 0.93
fyear2016 -10.67 0.88
fyear2017 -10.92 0.81
fyear2018 -10.00 0.88
fyear2019 -11.32 1.08
fyear2020 -10.27 0.76
fyear2021 -10.58 0.79
fyear2022 -11.35 0.89
fquarter4 -0.69 0.32
depth_sc -1.09 0.28
oxygen_sc 0.05 0.24
small_cod_density_sc 0.29 0.24
density_saduria_sc 0.43 0.26
flounder_density_sc -0.62 0.20
small_cod_density_sc:density_saduria_sc 0.09 0.20
density_saduria_sc:flounder_density_sc 0.07 0.14
Dispersion parameter: 0.13
Tweedie p: 1.47
Matérn range: 34.13
Spatial SD: 2.11
Spatiotemporal IID SD: 1.74
ML criterion at convergence: -731.195
See ?tidy.sdmTMB to extract these values as a data frame.
[1] "large cod"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably largeSpatiotemporal model fit by ML ['sdmTMB']
Formula: benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + small_cod_density_sc *
Formula: oxygen_sc + flounder_density_sc * oxygen_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
coef.est coef.se
fyear2015 -6.59 0.29
fyear2016 -6.48 0.22
fyear2017 -6.72 0.19
fyear2018 -6.56 0.22
fyear2019 -7.51 0.30
fyear2020 -6.71 0.18
fyear2021 -6.54 0.17
fyear2022 -6.79 0.21
fquarter4 0.83 0.12
depth_sc -0.24 0.07
small_cod_density_sc -0.09 0.06
oxygen_sc 0.21 0.07
flounder_density_sc 0.04 0.06
small_cod_density_sc:oxygen_sc 0.00 0.04
oxygen_sc:flounder_density_sc 0.04 0.04
Dispersion parameter: 0.27
Tweedie p: 1.63
Matérn range: 9.20
Spatiotemporal IID SD: 1.34
ML criterion at convergence: -7191.001
See ?tidy.sdmTMB to extract these values as a data frame.distinct: removed 2,968 rows (91%), 288 rows remaining
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.004954 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 49.54 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 401 [ 0%] (Warmup)
Chain 1: Iteration: 40 / 401 [ 9%] (Warmup)
Chain 1: Iteration: 80 / 401 [ 19%] (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%] (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%] (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%] (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%] (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%] (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%] (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%] (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%] (Warmup)
Chain 1: Iteration: 401 / 401 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 59.453 seconds (Warm-up)
Chain 1: 0.087 seconds (Sampling)
Chain 1: 59.54 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.004462 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 44.62 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 401 [ 0%] (Warmup)
Chain 1: Iteration: 40 / 401 [ 9%] (Warmup)
Chain 1: Iteration: 80 / 401 [ 19%] (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%] (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%] (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%] (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%] (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%] (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%] (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%] (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%] (Warmup)
Chain 1: Iteration: 401 / 401 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 28.885 seconds (Warm-up)
Chain 1: 0.04 seconds (Sampling)
Chain 1: 28.925 seconds (Total)
Chain 1: filter: removed 4 rows (80%), one row remainingmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'model' (character) with one unique value and 0% NAfilter: removed 3 rows (75%), one row remainingmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'model' (character) with one unique value and 0% NAmutate: new variable 'year' (double) with one unique value and 0% NA new variable 'fyear' (factor) with one unique value and 0% NA new variable 'fquarter' (factor) with one unique value and 0% NA new variable 'oxygen_sc' (double) with one unique value and 0% NA new variable 'depth_sc' (double) with one unique value and 0% NA new variable 'small_cod_density_sc' (double) with one unique value and 0% NAmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'xvar' (character) with one unique value and 0% NAmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'xvar' (character) with one unique value and 0% NAmutate: new variable 'year' (double) with one unique value and 0% NA new variable 'fyear' (factor) with one unique value and 0% NA new variable 'fquarter' (factor) with one unique value and 0% NA new variable 'oxygen_sc' (double) with one unique value and 0% NA new variable 'depth_sc' (double) with one unique value and 0% NA new variable 'flounder_density_sc' (double) with one unique value and 0% NAmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'xvar' (character) with one unique value and 0% NAmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'xvar' (character) with one unique value and 0% NAmutate: new variable 'species' (character) with one unique value and 0% NA new variable 'response' (character) with one unique value and 0% NA new variable 'sig' (character) with 2 unique values and 0% NAmutate: new variable 'species' (character) with one unique value and 0% NA new variable 'response' (character) with one unique value and 0% NA new variable 'sig' (character) with 2 unique values and 0% NAmutate: new variable 'group' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NAfilter: removed 7,107 rows (79%), 1,847 rows remaining[1] "small cod"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably largeSpatiotemporal model fit by ML ['sdmTMB']
Formula: saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc +
Formula: small_cod_density_sc * density_saduria_sc + flounder_density_sc *
Formula: density_saduria_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
coef.est coef.se
fyear2015 -8.73 1.10
fyear2016 -9.22 1.04
fyear2017 -9.89 0.99
fyear2018 -9.34 1.16
fyear2019 -8.25 1.15
fyear2020 -9.10 0.97
fyear2021 -9.70 0.96
fyear2022 -9.45 1.00
fquarter4 -1.47 0.36
depth_sc -0.98 0.30
oxygen_sc -0.40 0.25
small_cod_density_sc 0.44 0.39
density_saduria_sc 0.58 0.25
flounder_density_sc -0.56 0.20
small_cod_density_sc:density_saduria_sc 0.30 0.30
density_saduria_sc:flounder_density_sc 0.03 0.13
Dispersion parameter: 0.18
Tweedie p: 1.50
Matérn range: 65.62
Spatial SD: 2.40
Spatiotemporal IID SD: 1.27
ML criterion at convergence: -599.723
See ?tidy.sdmTMB to extract these values as a data frame.
[1] "small cod"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably largeSpatiotemporal model fit by ML ['sdmTMB']
Formula: benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + small_cod_density_sc *
Formula: oxygen_sc + flounder_density_sc * oxygen_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
coef.est coef.se
fyear2015 -6.01 0.31
fyear2016 -5.82 0.23
fyear2017 -5.78 0.19
fyear2018 -5.76 0.23
fyear2019 -6.34 0.31
fyear2020 -5.74 0.18
fyear2021 -5.98 0.17
fyear2022 -5.83 0.20
fquarter4 0.52 0.11
depth_sc -0.32 0.07
small_cod_density_sc 0.22 0.10
oxygen_sc 0.22 0.07
flounder_density_sc 0.03 0.07
small_cod_density_sc:oxygen_sc -0.16 0.08
oxygen_sc:flounder_density_sc 0.00 0.06
Dispersion parameter: 0.10
Tweedie p: 1.54
Matérn range: 14.30
Spatiotemporal IID SD: 0.92
ML criterion at convergence: -5234.112
See ?tidy.sdmTMB to extract these values as a data frame.distinct: removed 1,582 rows (86%), 265 rows remaining
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.005208 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 52.08 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 401 [ 0%] (Warmup)
Chain 1: Iteration: 40 / 401 [ 9%] (Warmup)
Chain 1: Iteration: 80 / 401 [ 19%] (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%] (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%] (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%] (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%] (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%] (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%] (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%] (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%] (Warmup)
Chain 1: Iteration: 401 / 401 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 214.791 seconds (Warm-up)
Chain 1: 0.771 seconds (Sampling)
Chain 1: 215.562 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.004676 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 46.76 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 401 [ 0%] (Warmup)
Chain 1: Iteration: 40 / 401 [ 9%] (Warmup)
Chain 1: Iteration: 80 / 401 [ 19%] (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%] (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%] (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%] (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%] (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%] (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%] (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%] (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%] (Warmup)
Chain 1: Iteration: 401 / 401 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 27.652 seconds (Warm-up)
Chain 1: 0.058 seconds (Sampling)
Chain 1: 27.71 seconds (Total)
Chain 1: filter: removed 4 rows (80%), one row remainingmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'model' (character) with one unique value and 0% NAfilter: removed 3 rows (75%), one row remainingmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'model' (character) with one unique value and 0% NAmutate: new variable 'year' (double) with one unique value and 0% NA new variable 'fyear' (factor) with one unique value and 0% NA new variable 'fquarter' (factor) with one unique value and 0% NA new variable 'oxygen_sc' (double) with one unique value and 0% NA new variable 'depth_sc' (double) with one unique value and 0% NA new variable 'small_cod_density_sc' (double) with one unique value and 0% NAmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'xvar' (character) with one unique value and 0% NAmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'xvar' (character) with one unique value and 0% NAmutate: new variable 'year' (double) with one unique value and 0% NA new variable 'fyear' (factor) with one unique value and 0% NA new variable 'fquarter' (factor) with one unique value and 0% NA new variable 'oxygen_sc' (double) with one unique value and 0% NA new variable 'depth_sc' (double) with one unique value and 0% NA new variable 'flounder_density_sc' (double) with one unique value and 0% NAmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'xvar' (character) with one unique value and 0% NAmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'xvar' (character) with one unique value and 0% NAmutate: new variable 'species' (character) with one unique value and 0% NA new variable 'response' (character) with one unique value and 0% NA new variable 'sig' (character) with 2 unique values and 0% NAmutate: new variable 'species' (character) with one unique value and 0% NA new variable 'response' (character) with one unique value and 0% NA new variable 'sig' (character) with 2 unique values and 0% NAmutate: new variable 'group' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NAfilter: removed 5,103 rows (57%), 3,851 rows remaining[1] "flounder"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably largeSpatiotemporal model fit by ML ['sdmTMB']
Formula: saduria_rel_weight ~ 0 + fyear + fquarter + depth_sc + oxygen_sc +
Formula: small_cod_density_sc * density_saduria_sc + flounder_density_sc *
Formula: density_saduria_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
coef.est coef.se
fyear2015 -8.68 1.14
fyear2016 -6.46 0.97
fyear2017 -8.57 0.97
fyear2018 -8.73 1.03
fyear2019 -8.42 1.07
fyear2020 -8.43 0.97
fyear2021 -11.38 1.27
fyear2022 -9.19 1.01
fquarter4 -0.57 0.21
depth_sc -0.33 0.21
oxygen_sc -0.04 0.22
small_cod_density_sc 0.41 0.13
density_saduria_sc 0.25 0.18
flounder_density_sc -0.20 0.15
small_cod_density_sc:density_saduria_sc 0.11 0.10
density_saduria_sc:flounder_density_sc -0.13 0.11
Dispersion parameter: 0.18
Tweedie p: 1.49
Matérn range: 74.50
Spatial SD: 2.20
Spatiotemporal IID SD: 1.59
ML criterion at convergence: -1790.092
See ?tidy.sdmTMB to extract these values as a data frame.
[1] "flounder"✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably largeSpatiotemporal model fit by ML ['sdmTMB']
Formula: benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc + small_cod_density_sc *
Formula: oxygen_sc + flounder_density_sc * oxygen_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
coef.est coef.se
fyear2015 -5.11 0.31
fyear2016 -4.94 0.24
fyear2017 -4.98 0.20
fyear2018 -5.50 0.24
fyear2019 -5.17 0.26
fyear2020 -4.70 0.19
fyear2021 -5.34 0.25
fyear2022 -5.03 0.21
fquarter4 0.04 0.10
depth_sc -0.36 0.10
small_cod_density_sc 0.22 0.06
oxygen_sc -0.05 0.07
flounder_density_sc -0.03 0.07
small_cod_density_sc:oxygen_sc 0.00 0.03
oxygen_sc:flounder_density_sc -0.02 0.05
Dispersion parameter: 0.17
Tweedie p: 1.50
Matérn range: 18.38
Spatial SD: 0.95
Spatiotemporal IID SD: 0.67
ML criterion at convergence: -5535.441
See ?tidy.sdmTMB to extract these values as a data frame.distinct: removed 3,590 rows (93%), 261 rows remaining
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.005659 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 56.59 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 401 [ 0%] (Warmup)
Chain 1: Iteration: 40 / 401 [ 9%] (Warmup)
Chain 1: Iteration: 80 / 401 [ 19%] (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%] (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%] (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%] (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%] (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%] (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%] (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%] (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%] (Warmup)
Chain 1: Iteration: 401 / 401 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 290.973 seconds (Warm-up)
Chain 1: 0.395 seconds (Sampling)
Chain 1: 291.368 seconds (Total)
Chain 1:
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.005409 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 54.09 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 401 [ 0%] (Warmup)
Chain 1: Iteration: 40 / 401 [ 9%] (Warmup)
Chain 1: Iteration: 80 / 401 [ 19%] (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%] (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%] (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%] (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%] (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%] (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%] (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%] (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%] (Warmup)
Chain 1: Iteration: 401 / 401 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 46.564 seconds (Warm-up)
Chain 1: 0.1 seconds (Sampling)
Chain 1: 46.664 seconds (Total)
Chain 1: filter: removed 4 rows (80%), one row remainingmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'model' (character) with one unique value and 0% NAfilter: removed 4 rows (80%), one row remainingmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'model' (character) with one unique value and 0% NAmutate: new variable 'year' (double) with one unique value and 0% NA new variable 'fyear' (factor) with one unique value and 0% NA new variable 'fquarter' (factor) with one unique value and 0% NA new variable 'oxygen_sc' (double) with one unique value and 0% NA new variable 'depth_sc' (double) with one unique value and 0% NA new variable 'small_cod_density_sc' (double) with one unique value and 0% NAmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'xvar' (character) with one unique value and 0% NAmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'xvar' (character) with one unique value and 0% NAmutate: new variable 'year' (double) with one unique value and 0% NA new variable 'fyear' (factor) with one unique value and 0% NA new variable 'fquarter' (factor) with one unique value and 0% NA new variable 'oxygen_sc' (double) with one unique value and 0% NA new variable 'depth_sc' (double) with one unique value and 0% NA new variable 'flounder_density_sc' (double) with one unique value and 0% NAmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'xvar' (character) with one unique value and 0% NAmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'xvar' (character) with one unique value and 0% NAmutate: new variable 'species' (character) with one unique value and 0% NA new variable 'response' (character) with one unique value and 0% NA new variable 'sig' (character) with 2 unique values and 0% NAmutate: new variable 'species' (character) with one unique value and 0% NA new variable 'response' (character) with one unique value and 0% NA new variable 'sig' (character) with 2 unique values and 0% NAmutate: new variable 'group' (character) with one unique value and 0% NA
mutate: new variable 'group' (character) with one unique value and 0% NANow do a separate model for adult cod looking at total prey
dd <- filter(d, group == "large cod")filter: removed 5,698 rows (64%), 3,256 rows remaining mesh <- make_mesh(dd,
xy_cols = c("X", "Y"),
cutoff = 5)
# Total model
# NOTE: turning off spatial here due to convergence and AIC
m_tot <- sdmTMB(tot_rel_weight ~ 0 + fyear + fquarter + depth_sc +
large_cod_density_sc*oxygen_sc + flounder_density_sc*oxygen_sc,
data = dd,
mesh = mesh,
family = tweedie(),
spatiotemporal = "IID",
spatial = "off",
time = "year")
sanity(m_tot)✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large print(m_tot)Spatiotemporal model fit by ML ['sdmTMB']
Formula: tot_rel_weight ~ 0 + fyear + fquarter + depth_sc + large_cod_density_sc *
Formula: oxygen_sc + flounder_density_sc * oxygen_sc
Mesh: mesh (isotropic covariance)
Time column: year
Data: dd
Family: tweedie(link = 'log')
coef.est coef.se
fyear2015 -4.15 0.21
fyear2016 -4.64 0.16
fyear2017 -4.53 0.13
fyear2018 -4.34 0.15
fyear2019 -4.76 0.22
fyear2020 -4.59 0.13
fyear2021 -4.36 0.13
fyear2022 -4.61 0.15
fquarter4 -0.19 0.10
depth_sc -0.01 0.05
large_cod_density_sc 0.16 0.06
oxygen_sc -0.03 0.05
flounder_density_sc -0.01 0.05
large_cod_density_sc:oxygen_sc 0.11 0.05
oxygen_sc:flounder_density_sc -0.09 0.04
Dispersion parameter: 0.60
Tweedie p: 1.71
Matérn range: 14.87
Spatiotemporal IID SD: 0.62
ML criterion at convergence: -7359.392
See ?tidy.sdmTMB to extract these values as a data frame. # Residuals
samps <- sdmTMBextra::predict_mle_mcmc(m_tot, mcmc_iter = 401, mcmc_warmup = 400)
SAMPLING FOR MODEL 'tmb_generic' NOW (CHAIN 1).
Chain 1:
Chain 1: Gradient evaluation took 0.004329 seconds
Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 43.29 seconds.
Chain 1: Adjust your expectations accordingly!
Chain 1:
Chain 1:
Chain 1: Iteration: 1 / 401 [ 0%] (Warmup)
Chain 1: Iteration: 40 / 401 [ 9%] (Warmup)
Chain 1: Iteration: 80 / 401 [ 19%] (Warmup)
Chain 1: Iteration: 120 / 401 [ 29%] (Warmup)
Chain 1: Iteration: 160 / 401 [ 39%] (Warmup)
Chain 1: Iteration: 200 / 401 [ 49%] (Warmup)
Chain 1: Iteration: 240 / 401 [ 59%] (Warmup)
Chain 1: Iteration: 280 / 401 [ 69%] (Warmup)
Chain 1: Iteration: 320 / 401 [ 79%] (Warmup)
Chain 1: Iteration: 360 / 401 [ 89%] (Warmup)
Chain 1: Iteration: 400 / 401 [ 99%] (Warmup)
Chain 1: Iteration: 401 / 401 [100%] (Sampling)
Chain 1:
Chain 1: Elapsed Time: 26.061 seconds (Warm-up)
Chain 1: 0.08 seconds (Sampling)
Chain 1: 26.141 seconds (Total)
Chain 1: mcmc_res <- residuals(m_tot, type = "mle-mcmc", mcmc_samples = samps)
dd$res <- as.vector(mcmc_res)
res_tot <- dd
# Range
range_tot <- tidy(m_tot, effects = "ran_pars") %>% filter(term == "range") %>% mutate(group = "large cod", model = "total")filter: removed 3 rows (75%), one row remainingmutate: new variable 'group' (character) with one unique value and 0% NA new variable 'model' (character) with one unique value and 0% NA # Spatial and spatiotemporal random effects
d_haul <- dd %>%
distinct(haul_id, .keep_all = TRUE)distinct: removed 2,968 rows (91%), 288 rows remaining preds_tot <- predict(m_tot, newdata = d_haul)
# Coefficients
coefs_tot <- bind_rows(tidy(m_tot, effects = "fixed", conf.int = TRUE)) %>%
mutate(species = "Cod (m)",
response = "Saduria",
sig = ifelse(estimate > 0 & conf.low > 0, "Y", "N"),
sig = ifelse(estimate < 0 & conf.high < 0, "Y", sig))mutate: new variable 'species' (character) with one unique value and 0% NA new variable 'response' (character) with one unique value and 0% NA new variable 'sig' (character) with 2 unique values and 0% NA coefs_tot <- coefs_tot %>% mutate(group = "large cod")mutate: new variable 'group' (character) with one unique value and 0% NA # Conditional effects: flounder
nd_flounder <- data.frame(expand_grid(
flounder_density_sc = seq(quantile(dd$flounder_density_sc, probs = 0.05),
quantile(dd$flounder_density_sc, probs = 0.95),
length.out = 50))) %>%
mutate(year = 2020,
fyear = as.factor(2020),
fquarter = as.factor(1),
oxygen_sc = 0,
depth_sc = 0,
large_cod_density_sc = 0)mutate: new variable 'year' (double) with one unique value and 0% NA new variable 'fyear' (factor) with one unique value and 0% NA new variable 'fquarter' (factor) with one unique value and 0% NA new variable 'oxygen_sc' (double) with one unique value and 0% NA new variable 'depth_sc' (double) with one unique value and 0% NA new variable 'large_cod_density_sc' (double) with one unique value and 0% NA preds_flounder_tot <- predict(m_tot, newdata = nd_flounder, re_form = NA, re_form_iid = NA, se_fit = TRUE)
pred_flounder_tot <- preds_flounder_tot %>% mutate(group = "large cod", xvar = "flounder")mutate: new variable 'group' (character) with one unique value and 0% NA new variable 'xvar' (character) with one unique value and 0% NAMake dataframes
coef_df <- bind_rows(bind_rows(coef_sad) %>% mutate(model = "Saduria"),
bind_rows(coef_ben) %>% mutate(model = "Benthos"))mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NAcoef_df <- coef_df %>% bind_rows(coefs_tot %>% mutate(model = "Total"))mutate: new variable 'model' (character) with one unique value and 0% NApred_cod_df <- bind_rows(bind_rows(pred_cod_sad) %>% mutate(model = "Saduria"),
bind_rows(pred_cod_ben) %>% mutate(model = "Benthos"))mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NApred_flounder_df <- bind_rows(bind_rows(pred_flounder_sad) %>% mutate(model = "Saduria"),
bind_rows(pred_flounder_ben) %>% mutate(model = "Benthos"))mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NAres_df <- bind_rows(bind_rows(res_sad) %>% mutate(model = "Saduria"),
bind_rows(res_ben) %>% mutate(model = "Benthos"))mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NAres_df <- res_df %>% bind_rows(res_tot %>% mutate(model = "Total"))mutate: new variable 'model' (character) with one unique value and 0% NArandom_df <- bind_rows(bind_rows(random_sad) %>% mutate(model = "Saduria"),
bind_rows(random_ben) %>% mutate(model = "Benthos"))mutate: new variable 'model' (character) with one unique value and 0% NA
mutate: new variable 'model' (character) with one unique value and 0% NArandom_df <- random_df %>% bind_rows(preds_tot %>% mutate(model = "Total"))mutate: new variable 'model' (character) with one unique value and 0% NArange_df <- bind_rows(range_tot, bind_rows(range_ben), bind_rows(range_sad))Plot spatial random effects
random_df <- random_df %>%
mutate(group = str_to_sentence(group))mutate: changed 1,916 values (100%) of 'group' (0 new NA)# Saduria
plot_map_fc +
geom_point(data = random_df %>% filter(model == "Saduria"),
aes(X*1000, Y*1000, color = omega_s), size = 0.9) +
scale_color_gradient2() +
facet_wrap(~factor(group, levels = c("Flounder", "Small cod", "Large cod"))) +
labs(color = "Spatial\nrandom effect") +
theme(axis.text.x = element_text(angle = 90),
legend.position = "bottom",
#legend.direction = "vertical",
legend.key.width = unit(0.6, "cm"),
legend.key.height = unit(0.2, "cm"))filter: removed 1,102 rows (58%), 814 rows remaining
ggsave(paste0(home, "/figures/supp/omega_sad.pdf"), width = 17, height = 9, units = "cm")
# Now do benthos (only for flounder)
plot_map_fc +
geom_point(data = random_df %>% filter(model == "Benthos" & group == "Flounder"),
aes(X*1000, Y*1000, color = omega_s), size = 0.9) +
scale_color_gradient2() +
facet_wrap(~ group) +
labs(color = "Spatial\nrandom effect") +
theme(axis.text.x = element_text(angle = 90),
legend.position = "right",
legend.direction = "vertical",
legend.key.width = unit(0.4, "cm"),
legend.key.height = unit(0.4, "cm"))filter: removed 1,655 rows (86%), 261 rows remaining
ggsave(paste0(home, "/figures/supp/omega_ben.pdf"), width = 11, height = 11, units = "cm")Plot spatiotemporal random effects
# Saduria
sad_eps_sc <- plot_map_fc +
geom_point(data = random_df %>% filter(model == "Saduria" & group == "Small cod"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
scale_color_gradient2() +
facet_wrap(~year) +
labs(color = "Spatiotemporal\nrandom effect") +
theme(legend.position = c(0.84, 0.16),
axis.text.x = element_text(angle = 90))filter: removed 1,651 rows (86%), 265 rows remainingWarning: A numeric `legend.position` argument in `theme()` was deprecated in ggplot2
3.5.0.
ℹ Please use the `legend.position.inside` argument of `theme()` instead.sad_eps_sc
ggsave(paste0(home, "/figures/supp/epsilon_sad_small_cod.pdf"), width = 17, height = 17, units = "cm")
sad_eps_lc <- plot_map_fc +
geom_point(data = random_df %>% filter(model == "Saduria" & group == "Large cod"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
scale_color_gradient2() +
facet_wrap(~ year) +
labs(color = "Spatiotemporal\nrandom effect") +
theme(legend.position = c(0.84, 0.16),
axis.text.x = element_text(angle = 90))filter: removed 1,628 rows (85%), 288 rows remainingsad_eps_lc
ggsave(paste0(home, "/figures/supp/epsilon_sad_large_cod.pdf"), width = 17, height = 17, units = "cm")
sad_eps_f <- plot_map_fc +
geom_point(data = random_df %>% filter(model == "Saduria" & group == "Flounder"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
scale_color_gradient2() +
facet_wrap(~ year) +
labs(color = "Spatiotemporal\nrandom effect") +
theme(legend.position = c(0.84, 0.16),
axis.text.x = element_text(angle = 90))filter: removed 1,655 rows (86%), 261 rows remainingsad_eps_f
ggsave(paste0(home, "/figures/supp/epsilon_sad_flounder.pdf"), width = 17, height = 17, units = "cm")
# Benthos
ben_eps_sc <- plot_map_fc +
geom_point(data = random_df %>% filter(model == "Benthos" & group == "Small cod"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
scale_color_gradient2() +
facet_wrap(~ year) +
labs(color = "Spatiotemporal\nrandom effect") +
theme(legend.position = c(0.84, 0.16),
axis.text.x = element_text(angle = 90))filter: removed 1,651 rows (86%), 265 rows remainingben_eps_sc
ggsave(paste0(home, "/figures/supp/epsilon_ben_small_cod.pdf"), width = 17, height = 17, units = "cm")
ben_eps_lc <- plot_map_fc +
geom_point(data = random_df %>% filter(model == "Benthos" & group == "Large cod"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
scale_color_gradient2() +
facet_wrap(~ year) +
labs(color = "Spatiotemporal\nrandom effect") +
theme(legend.position = c(0.84, 0.16),
axis.text.x = element_text(angle = 90))filter: removed 1,628 rows (85%), 288 rows remainingben_eps_lc
ggsave(paste0(home, "/figures/supp/epsilon_ben_large_cod.pdf"), width = 17, height = 17, units = "cm")
ben_eps_f <- plot_map_fc +
geom_point(data = random_df %>% filter(model == "Benthos" & group == "Flounder"), aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
scale_color_gradient2() +
facet_wrap(~ year) +
labs(color = "Spatiotemporal\nrandom effect") +
theme(legend.position = c(0.84, 0.16),
axis.text.x = element_text(angle = 90))filter: removed 1,655 rows (86%), 261 rows remainingben_eps_f
ggsave(paste0(home, "/figures/supp/epsilon_ben_flounder.pdf"), width = 17, height = 17, units = "cm")
# Total
tot_eps <- plot_map_fc +
geom_point(data = preds_tot, aes(X*1000, Y*1000, color = epsilon_st), size = 0.9) +
scale_color_gradient2() +
facet_wrap(~ year) +
labs(color = "Spatiotemporal\nrandom effect") +
theme(legend.position = c(0.84, 0.16),
axis.text.x = element_text(angle = 90))
tot_eps
Plot range
#pal <- brewer.pal(n = 8, name = "Dark2")[c(2, 7, 6)]
pal <- (brewer.pal(n = 11, name = "RdYlBu")[c(11, 4, 1)])
range_df %>% arrange(estimate)# A tibble: 7 × 5
term estimate std.error group model
<chr> <dbl> <dbl> <chr> <chr>
1 range 9.20 2.79 large cod benthos
2 range 14.3 4.14 small cod benthos
3 range 14.9 4.88 large cod total
4 range 18.4 4.36 flounder benthos
5 range 34.1 16.2 large cod saduria
6 range 65.6 30.1 small cod saduria
7 range 74.5 16.6 flounder saduriarange_df %>%
mutate(group = str_to_sentence(group),
model2 = ifelse(model == "benthos", "Benthic prey", model),
model2 = ifelse(model == "saduria", "Saduria", model2),
model2 = ifelse(model == "total", "Total prey", model2)) %>%
ggplot(aes(model2, estimate, color = factor(group, levels = c("Flounder", "Small cod", "Large cod")))) +
geom_point(size = 2) +
geom_hline(yintercept = 5, linetype = 2, alpha = 0.5) +
scale_color_manual(values = pal) +
labs(x = "", y = "Range (km)", color = "Predator") +
theme(aspect.ratio = 1,
legend.position = c(0.86, 0.86)) mutate: changed 7 values (100%) of 'group' (0 new NA) new variable 'model2' (character) with 3 unique values and 0% NA
ggsave(paste0(home, "/figures/supp/ranges.pdf"), width = 11, height = 11, units = "cm")Plot residuals
# Plot residuals
res_df %>%
mutate(group = str_to_title(group)) %>%
ggplot(aes(sample = res)) +
stat_qq(size = 0.75, shape = 21, fill = NA) +
facet_grid(model ~ group) +
stat_qq_line() +
labs(y = "Sample Quantiles", x = "Theoretical Quantiles") +
theme(aspect.ratio = 1)mutate: changed 21,164 values (100%) of 'group' (0 new NA)
ggsave(paste0(home, "/figures/supp/qq_relative_prey_weight.pdf"), width = 17, height = 17, units = "cm")Plot coefficients
coef_df$term %>% unique() [1] "fyear2015"
[2] "fyear2016"
[3] "fyear2017"
[4] "fyear2018"
[5] "fyear2019"
[6] "fyear2020"
[7] "fyear2021"
[8] "fyear2022"
[9] "fquarter4"
[10] "depth_sc"
[11] "oxygen_sc"
[12] "small_cod_density_sc"
[13] "density_saduria_sc"
[14] "flounder_density_sc"
[15] "small_cod_density_sc:density_saduria_sc"
[16] "density_saduria_sc:flounder_density_sc"
[17] "small_cod_density_sc:oxygen_sc"
[18] "oxygen_sc:flounder_density_sc"
[19] "large_cod_density_sc"
[20] "large_cod_density_sc:oxygen_sc" # Fix some names
coef_df2 <- coef_df %>%
filter(!grepl('year', term)) %>%
filter(!grepl('quarter', term)) %>%
mutate(term = str_remove_all(term, "_sc"),
term = str_remove_all(term, "density"),
term = str_replace_all(term, "_", ""),
term = str_replace_all(term, "geco", "ge co"),
term = str_replace_all(term, "llco", "ll co"),
term = str_replace(term, ":", " × "),
term = str_to_sentence(term),
group = str_to_sentence(group)) %>%
mutate(model2 = ifelse(model == "Saduria", "Prey=Saduria", NA),
model2 = ifelse(model == "Benthos", "Prey=Benthic prey", model2),
model2 = ifelse(model == "Total", "Prey=Total prey", model2)) %>%
mutate(sig2 = ifelse(sig == "Y", "N", "Y"))filter: removed 56 rows (52%), 52 rows remainingfilter: removed 7 rows (13%), 45 rows remainingmutate: changed 45 values (100%) of 'term' (0 new NA) changed 45 values (100%) of 'group' (0 new NA)mutate: new variable 'model2' (character) with 3 unique values and 0% NAmutate: new variable 'sig2' (character) with 2 unique values and 0% NAp1 <-
coef_df2 %>%
filter(model == "Benthos") %>%
ggplot(aes(estimate, term, alpha = sig2,
color = factor(group, levels = c("Flounder", "Small cod", "Large cod")))) +
geom_stripped_rows(aes(y = term), inherit.aes = FALSE) +
facet_wrap2(~model2, ncol = 2, scales = "free") +
geom_vline(xintercept = 0, linetype = 2, alpha = 0.5, color = "gray10", linewidth = 0.2) +
geom_point(position = position_dodge(width = 0.5), size = 1.5) +
geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0,
position = position_dodge(width = 0.5)) +
scale_alpha_manual(values = c(1, 0.4)) +
scale_color_manual(values = pal) +
labs(x = "", y = "", alpha = "95% CI crossing 0", color = "Predator") +
guides(color = guide_legend(title.position = "top", title.hjust = 0.5),
alpha = guide_legend(title.position = "top", title.hjust = 0.5)) +
theme(legend.position = c(0.75, 0.2),
legend.direction = "vertical",
legend.box.spacing = unit(-3, "pt"),
legend.margin = margin(0, 0, 0, 0)) +
NULLfilter: removed 27 rows (60%), 18 rows remainingp2 <-
coef_df2 %>%
filter(model == "Saduria") %>%
ggplot(aes(estimate, term, alpha = sig2,
color = factor(group, levels = c("Flounder", "Small cod", "Large cod")))) +
geom_stripped_rows(aes(y = term), inherit.aes = FALSE) +
facet_wrap2(~model2, ncol = 2, scales = "free") +
geom_vline(xintercept = 0, linetype = 2, alpha = 0.5, color = "gray10", linewidth = 0.2) +
geom_point(position = position_dodge(width = 0.5), size = 1.5) +
geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0,
position = position_dodge(width = 0.5)) +
scale_alpha_manual(values = c(1, 0.4)) +
scale_color_manual(values = pal) +
labs(x = "", y = "") +
guides(color = "none",
alpha = "none")filter: removed 24 rows (53%), 21 rows remainingp3 <-
coef_df2 %>%
filter(model == "Total") %>%
ggplot(aes(estimate, term, alpha = sig2,
color = "Large cod")) +
geom_stripped_rows(aes(y = term), inherit.aes = FALSE) +
facet_wrap2(~model2, ncol = 2, scales = "free") +
geom_vline(xintercept = 0, linetype = 2, alpha = 0.5, color = "gray10", linewidth = 0.2) +
geom_point(position = position_dodge(width = 0.5), size = 1.5) +
geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0,
position = position_dodge(width = 0.5)) +
scale_alpha_manual(values = c(1, 0.4)) +
scale_color_manual(values = pal[3]) +
labs(x = "", y = "") +
guides(color = "none",
alpha = "none")filter: removed 39 rows (87%), 6 rows remainingp2 / p1 / p3 +
plot_annotation(tag_levels = "a") +
plot_layout(axes = "collect", guides = "collect") &
theme(legend.position = "bottom") &
coord_cartesian(xlim = c(-1.55, 1.15))
ggsave(paste0(home, "/figures/coefs.pdf"), width = 11, height = 22, units = "cm")Plot year and quarter coefficients
# Fix some names
coef_df3 <- coef_df %>%
filter(grepl('year', term)) %>%
mutate(term = str_remove_all(term, "fyear"),
group = str_to_sentence(group),
term = as.numeric(term),
model2 = ifelse(model == "Benthos", "Benthic prey", model),
model2 = ifelse(model == "Saduria", "Saduria", model2),
model2 = ifelse(model == "Total", "Total prey", model2))filter: removed 52 rows (48%), 56 rows remainingmutate: converted 'term' from character to double (0 new NA) changed 56 values (100%) of 'group' (0 new NA) new variable 'model2' (character) with 3 unique values and 0% NAggplot(coef_df3, aes(term, exp(estimate), color = factor(group, levels = c("Flounder", "Small cod", "Large cod")),
fill = factor(group, levels = c("Flounder", "Small cod", "Large cod")))) +
facet_wrap(~model2, scales = "free", ncol = 1) +
#geom_line(position = position_dodge(width = 0.5)) +
geom_point(position = position_dodge(width = 0.4)) +
#geom_ribbon(aes(ymin = exp(conf.low), ymax = exp(conf.high)), alpha = 0.2, color = NA) +
geom_errorbar(aes(ymin = exp(conf.low), ymax = exp(conf.high)), alpha = 0.4, width = 0,
position = position_dodge(width = 0.4)) +
scale_color_manual(values = pal) +
scale_fill_manual(values = pal) +
labs(x = "Year", y = "Standardized coefficient", color = "") +
guides(color = guide_legend(title.position = "top", title.hjust = 0.5, ncol = 3),
fill = "none") +
theme(legend.position = c(0.5, 0.99),
legend.direction = "vertical",
legend.box.spacing = unit(-3, "pt"),
legend.margin = margin(0, 0, 0, 0),
strip.text.x.top = element_text(angle = 0, hjust = 0)) +
NULL
ggsave(paste0(home, "/figures/supp/coefs_year.pdf"), width = 11, height = 21, units = "cm")# Now do quarter
coef_df5 <- coef_df %>%
filter(term %in% c("fquarter4")) %>%
mutate(group = str_to_sentence(group),
model2 = ifelse(model == "Benthos", "Benthic prey", model),
model2 = ifelse(model == "Saduria", "Saduria", model2),
model2 = ifelse(model == "Total", "Total prey", model2)) %>%
mutate(sig2 = ifelse(sig == "Y", "N", "Y"))filter: removed 101 rows (94%), 7 rows remainingmutate: changed 7 values (100%) of 'group' (0 new NA) new variable 'model2' (character) with 3 unique values and 0% NAmutate: new variable 'sig2' (character) with 2 unique values and 0% NAggplot(coef_df5, aes(estimate, model2,
alpha = sig2,
color = factor(group, levels = c("Flounder", "Small cod", "Large cod")))) +
geom_vline(xintercept = 0, linetype = 2, alpha = 0.5, color = "gray10", linewidth = 0.2) +
geom_point(position = position_dodge(width = 0.5), size = 1.5) +
geom_errorbar(aes(xmin = conf.low, xmax = conf.high), width = 0,
position = position_dodge(width = 0.5)) +
scale_alpha_manual(values = c(1, 0.4)) +
scale_color_manual(values = pal) +
labs(x = "", y = "Quarter 4 effect", alpha = "95% CI crossing 0", color = "Predator") +
guides(color = guide_legend(title.position = "top", title.hjust = 0.5),
alpha = guide_legend(title.position = "top", title.hjust = 0.5)) +
geom_stripped_rows(aes(y = model2), inherit.aes = FALSE) +
theme(legend.position = "bottom",
legend.direction = "horizontal",
legend.box = "horizontal",
legend.box.spacing = unit(-3, "pt"),
legend.margin = margin(0, 0, 0, 0)) +
coord_cartesian(expand = 0)
ggsave(paste0(home, "/figures/supp/coefs_quarter.pdf"), width = 17, height = 11, units = "cm")Conditional effects
# Which CI?
# https://www.calculator.net/confidence-interval-calculator.html
pred_df <- bind_rows(pred_cod_df, pred_flounder_df) %>%
mutate(group = str_to_sentence(group),
sad = ifelse(density_saduria_sc == min(density_saduria_sc), "Low", NA),
sad = ifelse(density_saduria_sc == max(density_saduria_sc), "High", sad)) %>%
drop_na(sad)mutate: changed 1,500 values (100%) of 'group' (0 new NA) new variable 'sad' (character) with 3 unique values and 20% NAdrop_na: removed 300 rows (20%), 1,200 rows remainingpred_df2 <- bind_rows(pred_flounder_df,
pred_flounder_tot %>% mutate(model = "Total")) %>%
mutate(group = str_to_sentence(group),
density_saduria_sc = replace_na(density_saduria_sc,
median(density_saduria_sc, na.rm = TRUE)))mutate: new variable 'model' (character) with one unique value and 0% NAmutate: changed 50 values (5%) of 'density_saduria_sc' (50 fewer NA) changed 950 values (100%) of 'group' (0 new NA)# 75% CI!!
ggplot(pred_df %>% filter(model == "Saduria" & xvar == "flounder"),# & !group == "Large cod"),
aes(flounder_density_sc, exp(est), color = sad, fill = sad)) +
geom_ribbon(aes(ymin = exp(est - 1.150*est_se), ymax = exp(est + 1.150*est_se)),
alpha = 0.3, color = NA) +
geom_line() +
facet_wrap(~factor(group, levels = c("Flounder", "Small cod", "Large cod")),
scales = "free",
ncol = 3) +
scale_color_brewer(palette = "Dark2") +
scale_fill_brewer(palette = "Dark2") +
coord_cartesian(#xlim = c(-0.5, 0.5),
expand = 0, clip = "off") +
labs(x = "Flounder density", y = "Relative saduria weight",
color = "Saduria", fill = "Saduria") +
theme(legend.position = c(0.95, 0.84),
strip.text.x.top = element_text(angle = 0, hjust = 0)) +
scale_x_continuous(breaks = c(-1, 0, 1)) +
NULLfilter: removed 900 rows (75%), 300 rows remaining
ggsave(paste0(home, "/figures/conditional_saduria_flounder.pdf"), width = 19, height = 7, units = "cm")Showing conditional effects of oxygen on small cod feeding on benthos
dd <- filter(d, group == "small cod")filter: removed 7,107 rows (79%), 1,847 rows remaining mesh <- make_mesh(dd,
xy_cols = c("X", "Y"),
cutoff = 5)
# Benthic model
m_ben <- sdmTMB(benthic_rel_weight ~ 0 + fyear + fquarter + depth_sc +
small_cod_density_sc*oxygen_sc + flounder_density_sc*oxygen_sc,
data = dd,
mesh = mesh,
family = tweedie(),
spatiotemporal = "IID",
spatial = "off",
time = "year")
sanity(m_ben)✔ Non-linear minimizer suggests successful convergence✔ Hessian matrix is positive definite✔ No extreme or very small eigenvalues detected✔ No gradients with respect to fixed effects are >= 0.001✔ No fixed-effect standard errors are NA✔ No standard errors look unreasonably large✔ No sigma parameters are < 0.01✔ No sigma parameters are > 100✔ Range parameter doesn't look unreasonably large nd <- data.frame(oxygen = seq(quantile(d$oxygen, probs = 0.05), quantile(d$oxygen, probs = 0.95),
length.out = 50)) %>%
mutate(year = 2020,
fyear = as.factor(2020),
fquarter = as.factor(1),
density_saduria_sc = 0,
flounder_density_sc = 0,
depth_sc = 0,
small_cod_density_sc = 0) %>%
mutate(oxygen_sc = (oxygen - mean(d$oxygen)) / sd(d$oxygen))mutate: new variable 'year' (double) with one unique value and 0% NA new variable 'fyear' (factor) with one unique value and 0% NA new variable 'fquarter' (factor) with one unique value and 0% NA new variable 'density_saduria_sc' (double) with one unique value and 0% NA new variable 'flounder_density_sc' (double) with one unique value and 0% NA new variable 'depth_sc' (double) with one unique value and 0% NA new variable 'small_cod_density_sc' (double) with one unique value and 0% NAmutate: new variable 'oxygen_sc' (double) with 50 unique values and 0% NA p <- predict(m_ben, newdata = nd, re_form = NA, re_form_iid = NA, se_fit = TRUE)
ggplot(p, aes(oxygen, exp(est))) +
geom_line() +
theme_sleek(base_size = 14) +
geom_hline(yintercept = 0.0025, col = "red") +
geom_hline(yintercept = 0.0040, col = "red") +
geom_vline(xintercept = 4.8, col = "red") +
geom_vline(xintercept = 7.6, col = "red") +
NULL